Impact of different destocking strategies on the resilience of dry rangelands

Abstract Half of the world's livestock live in (semi‐)arid regions, where a large proportion of people rely on animal husbandry for their survival. However, overgrazing can lead to land degradation and subsequent socio‐economic crises. Sustainable management of dry rangeland requires suitable stocking strategies and has been the subject of intense debate in the last decades. Our goal is to understand how variations in stocking strategies affect the resilience of dry rangelands. We describe rangeland dynamics through a simple mathematical model consisting of a system of coupled differential equations. In our model, livestock density is limited only by forage availability, which is itself limited by water availability. We model processes typical of dryland vegetation as a strong Allee effect, leading to bistability between a vegetated and a degraded state, even in the absence of herbivores. We study analytically the impact of varying the stocking density and the destocking adaptivity on the resilience of the system to the effects of drought. By using dynamical systems theory, we look at how different measures of resilience are affected by variations in destocking strategies. We find that the following: (1) Increasing stocking density decreases resilience, giving rise to an expected trade‐off between productivity and resilience. (2) There exists a maximal sustainable livestock density above which the system can only be degraded. This carrying capacity is common to all strategies. (3) Higher adaptivity of the destocking rate to available forage makes the system more resilient: the more adaptive a system is, the bigger the losses of vegetation it can recover from, without affecting the long‐term level of productivity. The first two results emphasize the need for suitable dry rangeland management strategies, to prevent degradation resulting from the conflict between profitability and sustainability. The third point offers a theoretical suggestion for such a strategy.


| INTRODUC TI ON
Dry rangelands are arid, semi-arid, and dry subhumid ecosystems that are dedicated mainly to livestock and wildlife grazing and browsing. They cover roughly a third of the Earth's land surface (ILRI et al., 2021;Safriel et al., 2006) and are home to about half of the world's livestock (Safriel et al., 2006). Water scarcity, due to a low-precipitation-to-evapotranspiration ratio, limits the possibility of growing crops (Hobbs et al., 2008). Still, dry rangelands support the livelihoods of hundreds of millions of the world's poorest people through animal husbandry (de Haan, 2016;Herrero et al., 2009;Naess & Bardsen, 2013). This pivotal role in global food security is expected to increase, due to growth in both drylands' human population and their per capita consumption of meat and other animal products (de Haan, 2016;Herrero et al., 2015). Furthermore, dry rangelands provide ecosystem services such as carbon storage, soil formation, flooding control, and maintenance of biodiversity, as well as aesthetic and cultural value (FAO, 2009;Henderson et al., 2015;Sala et al., 2017;Sandhage-Hofmann, 2016).
The first hazard dry rangelands regularly experience, independently of how they are managed, is drought. This extreme water shortage has a direct negative impact on the health and growth of dryland vegetation (Gouveia et al., 2017;NOAA et al., 2022;Vicente-Serrano et al., 2013;Wang et al., 2015), which can be exacerbated by herbivory (Fritts et al., 2018), with subsequent negative impact on domestic herds (Catley et al., 2013;Catley et al., 2014;Naess & Bardsen, 2013). Another threat dry rangelands experience is overgrazing, which can lead to soil and vegetation degradation (Gonzalez & Ghermandi, 2021). In such contexts, drought can trigger regular collapses in animal numbers (Coppock et al., 2008;Desta & Coppock, 2002). Grazing pressure is currently intensifying due to demographic growth and other social factors affecting drylands, such as changing land use and tenure, sedentarisation of mobile pastoralists, changes in herd size and structure, and increased use of supplementary feed (Reid et al., 2014;Thornton et al., 2009). Overgrazing, in combination with drought, has been identified as a cause of 'desertification' (Brandt & Thornes, 1996;Geist & Lambin, 2004;Vetter, 2009;Yassoglou et al., 2017), a hardly reversible loss of productivity of dry rangelands that is also known as land degradation. There is a high confidence in the fact that desertification and climate change will cause future reductions in crop and livestock productivity (Mirzabaev et al., 2019). Accordingly, the prevention of dry rangeland degradation is a great matter of concern and research (Briske et al., 2020;Campbell et al., 2006;Jakoby et al., 2015;Sandford & Scoones, 2006;Tietjen & Jeltsch, 2007;Vetter, 2005).
A key notion in the study of dryland degradation is that of resilience of an ecosystem, that is, its capacity to maintain itself in a desirable ecological state when subjected to perturbations (Holling, 1973). In socio-ecological systems such as dry rangelands, it is important to understand how human decisions impact the system's capacity to recover from uncontrolled perturbations.
In particular, the expected increase in the frequency, duration, and intensity of droughts in many parts of the world (Cherlet et al., 2018;Vetter, 2009) calls for a better understanding of the strategies that enhance dry rangelands' resilience to drought.
In practice, different preparation and response strategies to drought are commonly implemented, with the two main strategies (Torell et al., 2010) consisting of either maintaining low stocking densities or having an adaptive (also called 'flexible', 'opportunistic', or 'tracking') stocking density. The former aims at keeping the number of animals low to minimize the risk of overstocking during dry events.
The latter involves actively adjusting the herd size according to rainfall or available forage. Choosing whether to apply a low-density or an adaptive stocking strategy in dry rangelands has been at the core of animated debates for decades (Campbell et al., 2006;Sandford & Scoones, 2006;Torell et al., 2010), which have not been settled yet (Sandhage-Hofmann, 2016). These debates have given rise to a more theoretical argument on whether dry rangelands follow equilibrium or non-equilibrium dynamics (see Briske et al., 2020;Vetter, 2005 for an overview of this issue). It is worth mentioning that the lowdensity and adaptive strategies have erroneously been presented as the two opposing extremes of a spectrum, when in reality, intensity and adaptivity of stocking density are two separate axes of variation in livestock management strategies (Campbell et al., 2006), as shown in Figure 1.
This means that changes in density and adaptivity not only have distinct implications in terms of productivity and resilience of the rangeland but also that they can be combined. Another shortcoming identified by some of the actors in the debate is the lack of generality of previous studies, with many results being extrapolated from case studies (Campbell et al., 2006) or tied to particular conditions (Sandford & Scoones, 2006). Addressing these gaps, this work investigates the impact of different combinations of stocking density and adaptivity on the resilience of a generic dry rangeland using mathematical modeling.

F I G U R E 1
The four extreme management strategies and how they relate to our control parameters and . Past research on optimal grazing strategies in dry rangelands typically opposed the poles II and IV, where II was associated with traditional pastoral systems and IV with commercial pasture management (Campbell et al., 2006). In this work, decreasing the parameter increases the stocking density and increasing increases the stocking adaptivity.
Our model is similar to that of van de Koppel and Rietkerk (2000), which is a coupled system of ordinary differential equations (ODEs) based on pioneering work on the stability of consumer-resource systems (Noy-Meir, 1975;Rosenzweig & MacArthur, 1963). In contrast with these earlier models, by taking into account dryland vegetation's vulnerability to land degradation, the van de Koppel model predicts the possibility of an almost irreversible loss of vegetation and the subsequent collapse of herbivores in the system. In analogous models for temperate rangelands, such an abrupt irreversible change can be observed only if limitations other than forage availability are put on herbivore growth (Noy-Meir, 1975). Van de Koppel et al. show that these collapses to a degraded state are less likely to happen in systems where herbivores and vegetation are coupled than in systems where the two biomasses are uncoupled (van de Koppel & Rietkerk, 2000). They conclude that adaptive management strategies may not preclude the irreversible collapse of the system if they are not rapid enough to prevent soil degradation. In a spatial simulation-based study designed for semi-arid rangelands, Jakoby et al. (2015) compare two stocking strategies in which herbivore growth depends fully on vegetation intake, with a fixed stocking goal in one case and a goal that depends on available forage in the other case. Their simulations suggest that a system with a high constant stocking goal is more prone to collapse than a system where the stocking goal is high but adapts to the available forage. However, even though both van de Koppel and Rietkerk (2000) and Jakoby et al. (2015) show how adaptivity of the stocking density can help avoid catastrophic shifts, neither investigates the resilience of the system to external perturbations, such as droughts.
In Fletcher and Hilbert (2007), a generic consumer-resource model applicable to temperate rangelands is analyzed numerically to assess how different management choices can affect resilience to external perturbations. The authors show that different management strategies can yield the same long-term level of productivity but drastically differ in terms of the resilience of the system. Like van de Koppel and Rietkerk (2000), Fletcher and Hilbert find that totally decoupling herbivore growth from vegetation can be detrimental to the vegetation and, hence, to the whole system. However, the strategies compared by Fletcher and Hilbert (2007) aggregate several types of undefined management actions with no clear mapping to real-world processes. Moreover, animal dynamics are not modeled explicitly and do not incorporate the impact of forage availability on herbivore growth. In contrast, the present work targets a single type of managerial action that regulates herd size. In our model, herbivore growth depends exclusively on the amount of vegetation available (i.e., we do not consider supplementary feed). Therefore, the herd size is actively managed only through the removal of animals or their (re)introduction into the system. Our model focuses exclusively on destocking, the reduction of the number of animals in the system.
In practice, destocking is being increasingly used to mitigate herbivory pressure on dryland vegetation, minimize livestock dieoffs during droughts, and generate positive socio-economic effects (Morton & Barton, 2002). A great proportion of ranchers use destocking as one of their strategies to respond to droughts (Kachergis et al., 2014;Salmoral et al., 2020), while humanitarian programs in dry, famine-prone areas regularly facilitate emergency sales and slaughtering of livestock during droughts (Abebe et al., 2008;Aklilu & Wekesa, 2002;Morton & Barton, 2002). Existing destocking strategies differ mainly in terms of their baseline rate and their degree of adaptivity to environmental changes (Morton & Barton, 2002). As the baseline destocking rate determines the herd size, comparing these strategies links us back to the debate around low-density versus adaptive strategies. Our aim is to investigate how both of these strategies affect the resilience of dry rangeland systems.
To our knowledge, no general model exists that explicitly defines different destocking strategies and compares their impact on dry rangelands. In this work, we formulate a broadly applicable mathematical model and explain how variations in the destocking rate impact the productivity and the resilience of the desirable steady state of a generic dry rangeland system subject to drought. Avoiding the classical polarization of the debate, we look separately at the effects of variations in stocking density and stocking adaptivity. In terms of productivity, we focus on the long-term herd size and how it varies with changes in management parameters. In terms of resilience, we investigate how different destocking strategies affect the existence of a sustainable productive state and the threshold before a drought-induced catastrophic collapse. We consider the possibility of, and the potential mechanisms for increasing the resilience of a dry rangeland without reducing its long-term productivity.

| Model
Similarly to the model in van de Koppel and Rietkerk (2000) and the general predator-prey model in , we represent the plants and herbivores dynamics with a comprehensive consumerresource system of two coupled ODEs (1) with generic functional forms under simple constraints (a1-a6, detailed later). Conducting our analysis on such a general form means that our results hold for all the particular cases it covers. Our findings are, hence, more robust and general than the ones derived from a specific formulation.
The non-negative state variables V and H are the vegetation and herbivore densities, respectively. Note that H is usually called 'stocking rate' in rangeland management. f(V) is the vegetation growth rate. g(V) is the per capita consumption rate of vegetation by animals. is the plant-to-animal conversion factor, while D ; (V) is the per capita animal loss rate that specifies how fast animals are removed from the system and how this removal varies with and ≥ 0. Typically, in consumer-resource (predator-prey) models, this loss term is a constant. In our model, is the destocking adaptivity parameter that represents the degree of coupling between the animal loss rate and the density of forage available.
When = 0, the loss term D ; (V) reduces to a constant, whose value is determined by . Hence, fixes the baseline destocking rate that sets the intensity of the per capita loss rate D ; (V) and by doing so, fixes the long-term stocking density. In cases where > 0, the parameter still fixes the baseline destocking rate and the long-term stocking density but D ; (V) will increase as vegetation density declines. The degree of this adaptation, consisting of removing more animals as the available forage decreases, depends on : the greater the value of , the greater the adaptivity of D ; (V) to lack of forage. In practice, an increased animal mortality rate correlated to forage scarcity can be due to either natural causes (such as dehydration, illnesses, lower immunity, etc.) or a management strategy. Even though the term D ; (V) can be interpreted as a constant or resource-dependent loss of animals due to natural causes, we refer to it as a 'destocking strategy' in the rest of the text. We give a more precise form for D ; (V) in constraint a4. Table 1 compiles the constraints we are adding to the system of equations (1). The set of constraints a1 applies to the plant growth function f(V). We assume that there is a strong feedback, typical of water scarce environments, between vegetation and soil water, leading to a strong Allee effect; that is, the existence of a critical density threshold below which the population growth is negative.

| Constraints to the functional forms of f, g and D
This is not the case for all dry rangelands, as discussed in van de Koppel and Rietkerk (2000). Ecologically, the strong Allee effect in drylands corresponds to the combined phenomena of plantplant facilitation (Kéfi et al., 2007;Rietkerk et al., 2004) and selfreinforcement of the bare ground (Saco et al., 2007): in (semi-) arid regions, existing plants facilitate the growth and survival of those nearby by locally reducing evaporation and enhancing water infiltration (Davies et al., 2007;Holmgren et al., 1997;Holzapfel et al., 2006) (see Callaway (2007) for a description of the main plant-plant facilitation mechanisms), whereas bare soils are subject to increased erosion, which prevents the establishment of new plants (Saco et al., 2007). These feedback loops consolidate the bare ground state in such a way that plant growth and establishment are difficult at low vegetation densities (Courchamp, 2008).
As mentioned earlier, in the case of a strong Allee effect, vegetation growth is negative below a critical density (as for example in Figure 2), called the Allee threshold. Similar to the vegetation carrying capacity K, the Allee threshold A is typically a composite parameter that depends on various environmental conditions as well as the plant species. We will see in the Results (Section 3) that the existence of this positive feedback gives rise to bistability in vegetation productivity and is, hence, associated with the possibility of a discontinuous transition between alternative steady states, even in the absence of herbivores. An abrupt transition between a productive and a degraded state corresponds to the long-term degradation of the rangeland. Figure 2 gives an example of a function f(V) that satisfies a1.
Constraints a2 and a3 on the per capita consumption rate g of the herbivores ensure that the animals feed on existing vegetation and that a higher density of forage makes consumption more efficient: as the animals do not have to move as much to find food, their foraging efficiency is improved. It is noteworthy that the widely used linear, Holling II, and Holling III functional responses satisfy requirements a2 and a3 on g(V).
Constraint a4 and its subconstraints a5 and a6 describe the management strategies we are studying. Our goal is to understand the influence on resilience of varying destocking strategies across the two axes of variation presented in Figure 1. The constraints on the destocking rate D ; allow us to investigate separately the effects of the stocking density at equilibrium and of the adaptivity to changes in vegetation density. The control parameter sets the stocking density. By constraint a5, the destocking rate increases as increases.
Therefore, as can be seen in the results (Section 3), the stocking density decreases as increases. We will also see that the vegetation density at equilibrium V eq depends on . The expression Vegetation growth is continuous, A is its Allee threshold, K is its carrying capacity. a2 There is no foraging in the absence of vegetation.
sets the destocking adaptivity such that for a fixed , all yield the same equilibrium. a5 X( ) > 0, X � ( ) ≥ 0 sets the stocking density. There is no restocking: we consider only naturally rebuilding herds.

a6
< < 1 is a negligible term that prevents division by 0.
compares the vegetation density at equilibrium V eq to the current vegetation density V: the lower V is, the greater Note that we assume that the per capita animal loss rate D ; (V) is purely a management decision and that we neglect the loss of animal biomass due to metabolic expenses, to simplify our analysis. Importantly, in the case of a constant metabolic expense rate, as is generally found in the literature, adding such a term, that is, rewriting D ; (V) as m + D ; (V), with m > 0 does not affect our results, provided that that is, provided that the animal growth rate when vegetation is unlimited is greater than the metabolic expense rate.
In addition to constraints a1-a6, there are several implicit constraints that are already enclosed in the system of coupled equations (1): • Because animal growth depends on vegetation only, there is no other source of calories for the livestock, that is, no supplementary feed.
• Herd growth is limited by food availability only, not by space limitation or other density-dependent factors. This assumption is motivated by the fact that dry rangelands are typically extensive exploitation systems, where herd growth is limited by forage availability, which is itself limited by water availability.
• The dynamics described are spatially homogeneous and continuous in time, with constant environmental parameters A and K . The stable steady states represent the possible long-term configurations of the system. Perturbations of the state variables can be applied to represent the effects of below and above average conditions, as well as exceptional events that are external to the model. This is discussed in more detail in the 'resilience' paragraph below. The model can apply to ranchers and pastoralists but ignores herd mobility by averaging the dynamics over space. Because it also assumes uniform management, that is, a single given destocking strategy over the whole domain, the model is suitable for one management unit, such as one ranch or one uniformly managed pastoral area. The model can, hence, apply to a wide range of unit sizes, from tens to thousands of hectares.

| Analysis
Our goal is to compare the resilience of a productive system under different degrees of stocking density and adaptivity, that is, for different values of the stocking density parameter and of the adaptivity parameter . For this, we first need to understand under which conditions a system is productive and how to measure its resilience to the effects of drought, which we model as a perturbation as described in the following. Then, we can compare two systems with identical parameters except for the control parameter of interest, namely or . Similarly to classical work on consumer-resource systems (Noy- Meir, 1975;Rosenzweig & MacArthur, 1963), we use linear stability and graphical analysis of the phase plane (defined below) to both understand the range of possible behaviors of the system and define resilience metrics.

| Phase plane analysis
The phase plane is the space of all possible states (V; H) of the system and their evolution, for a given parametrisation. It features the solution trajectories and, in particular, the steady states of the system (1). Each point (V; H) can represent an initial condition of (1) that will evolve following a unique solution trajectory governed by the equations. We use the phase plane to summarize the stability results and explain: • in which cases vegetation and herbivory are theoretically incompatible and why this is so.
• how different destocking strategies affect the resilience of a productive state.
To draw the phase plane, we need to consider the nullclines, that is, the sets of points for which one of the two populations does not change. These curves, solutions to the equations dV dt = 0 and dH dt = 0, (which yields the two vegetation and the two herbivore nullclines, respectively), partition the phase plane into areas of different behaviors. The steady states of the system lie at the intersections of a vegetation and herbivore nullcline.
Example of vegetation growth function satisfying the constraint a1. The function given here is of the shape where r is a positive constant.

| Sustainable productive state
A sustainable productive state of the system is defined as a steady state V eq ; H eq that is stable and strictly positive in its two components. We refer to its H-component, H eq , as the 'long-term stocking density' or the 'long-term stocking goal' of the system. Although positive stable limit cycles are possible, we exclude them from our analysis, as managers are unlikely to try to achieve a fluctuating herd size.

| Resilience of a sustainable productive state
When the system is at a sustainable productive state, it can either recover fully or shift to a less desirable equilibrium after a state variable perturbation. As mentioned previously, there exists a threshold in the phase plane, beyond which the system cannot recover to the stable productive state. This threshold, also known as the 'separatrix', is the boundary between the degraded and the productive states' basins of attraction (see Figure 3). In other words, the long-term impact of a perturbation on the system depends on whether or not it has sent the state variables outside the basin of attraction of the stable productive state (van de Koppel & Rietkerk, 2000;van Voorn et al., 2007). The position of the separatrix depends on the values of the parameters and, hence, on the management strategy adopted. In extreme cases, a change in parameters can cause a stable steady state to lose its stability. This critical event corresponds to a bifurcation of the system. Importantly, the boundaries of a stable state's basin of attraction can be used to design metrics of that state's resilience to perturbations (Dakos & Kéfi, 2022;Krakovská et al., 2021), as explained in the next paragraph.

| Resilience metrics
The concept of resilience, although central in dynamical systems, is not unequivocally characterized: there exist a variety of definitions and measures of a system's resilience (Dakos & Kéfi, 2022;Krakovská et al., 2021;Rotz & Fraser, 2015). We define here two different measures of the resilience of the desirable steady state for given management parameters and . Each of the measures captures a different aspect of resilience.
• The distance to bifurcation is a measure of the resilience of the system to changes in parameter values. It is defined as the minimum amount of change needed along a given parameter to cause the disappearance of the desirable (sustainably productive) steady state (Dakos & Kéfi, 2022). It, therefore, measures how much management strategies are allowed to change, even gradually, before the system would collapse. We study it for both parameters. First, for changes in the parameter , we write dist bif ( ) ,

defined, for a fixed value, as
where is the current value and bif is the value of at which the closest bifurcation occurs. When considering changes in the parameter , we write dist bif ( ). It is defined for a fixed , as where is the current value and bif is the value of at which the closest bifurcation occurs.
• The maximal loss ratio is a measure of the resilience of the desirable steady state to external perturbations. It reflects how many times greater the vegetation density at equilibrium is, Qualitative 3D stability landscape of our two dimensional system. The separatrix (green continuous line) marks the limit between the basin of attraction of the sustainable productive steady state V eq ; H eq and the basin of attraction of the degraded state (0;0). The position of the separatrix depends on the destocking strategy.
compared with the vegetation density at which the tipping point to the degraded state occurs. V pmax ( ; ) is the minimal value of V still in the basin of attraction of the productive steady state, when no change in H is applied (see below for a description of the perturbation and Figure 3 for an illustration of V pmax ). In practice, V pmax ( ; ) is the lowest the vegetation can get, due to an external perturbation (such as a drought), without the system collapsing to the degraded state.

| Modeling drought as a perturbation
We model the effects of drought as a sudden reduction of vegetation, with no change in H. This assumes a rapid effect of drought on vegetation (Noy-Meir, 1975;Zhao et al., 2020). The perturbed state, therefore, has the form V pert ; H eq , where the perturbed vegeta- The perturbation is impulsive (pulse perturbation) in the sense that we set the initial conditions of the system (1) to a perturbed state, then the system is subject to its usual dynamics. The perturbation is isolated in the sense that we always allow the system to recover or degrade fully after the perturbation is applied.

| RE SULTS
In this section, we first describe the range of behaviors

| General behavior of the system (phase plane analysis)
Our model is a general consumer-resource (predator-prey) system with a strong Allee effect on the resource. Equations of this type have been studied analytically in  for the case where the consumer's mortality rate (i.e., our destocking rate) is constant, that is, = 0. Our analysis showed that the behavior of the . The nullcline ℋ(V) splits the phase plane into the area above it, where vegetation density decreases, and the area below it, where vegetation density increases.
Finally, we have the non-trivial herbivore nullcline D ; This equation is independent of H because our model assumes that growth and destocking both depend linearly on herbivore density. From the shapes of the non-trivial nullclines, we see that, for a given parametrisation, there can be at most one sustainable productive steady state, V eq ; H eq . This stable node or sink is a potential management goal, with H eq the long-term stocking target of the rangeland manager. Note that with our definition of D ; , for any given , the management nullclines coincide for all values . On the other hand, varying the value of shifts the position of the nullcline: toward the right as we increase and toward the left as is decreased.

Moreover
Hence, the values of V eq and H eq as well as the stability properties of the four steady states (0; 0), (A; 0), (K; 0), and V eq ; H eq depend on the value of , while they remain unchanged by variations in . There are five broad categories of behavior: region (1) V eq > K. Because it implies that H eq is negative and that the equilibrium V eq ; H eq is unstable, this configuration is not ecologically relevant. Depending on the initial conditions: either the system settles into to the vegetation only equilibrium (K; 0), or it degrades itself to (0; 0) if the initial density of herbivores was too high and/or the initial plant density too low. region (2) V crit < V eq < K, where V crit is the maximal argument of ℋ(V), that is, ℋ V crit is the maximum value of the function ℋ(V) .
Therefore, ℋ V crit is by definition the animal carrying capacity H max , that is, the maximal herbivore density such that the productive steady state is stable. Importantly, H max is a fixed value that does not depend on the management function D ; (V) but only on the biological functions f(V) and g (V). Because the equilibrium V eq ; H eq is strictly positive and stable, this is a desirable configuration. The system will stabilize at the productive state V eq ; H eq , provided that the initial conditions lie within its basin of attraction. As mentioned earlier, the animal density at equilibrium H eq can be thought of as a long-term 'stocking goal' of the rangeland manager. The higher the stocking goal, the lower the vegetation density at equilibrium V eq , until the threshold value V crit ; H max is met. This point is a Hopf bifurcation, meaning a point of critical transition in the behavior of the system, where the productive equilibrium loses its stability and (stable or unstable) limit cycles appear. Such bifurcation at the top of the vegetation nullcline hump is not specific of a system with a strong Allee effect. It has been shown for other consumer-resource systems (e.g., Hilker & Schmitz, 2008).

region (3) Distinguishing region 3 from region 4 is relevant only
if the Hopf bifurcation abovementioned is supercritical, that is, gives rise to stable limit cycles. When the Hopf bifurcation is instead subcritical, region 3 does not exist and the dynamics are directly the ones described for region 4. Whether or not the bifurcation is supercritical or subcritical and, therefore, whether or not the oscillations are stable or unstable, depends on the specification of g and f. The details of the conditions for the Hopf cycle to be supercritical and, hence, for the existence of a region 3 are given in . When this is the case, Immediately after the Hopf bifurcation, for each value of , there is a range V ; V crit of values of V eq such that the system admits a stable limit cycle. Even though the productive steady state V eq ; H eq is now unstable, the system is still productive, with vegetation and animal densities oscillating in time, provided that the initial values lie in the basin of attraction of this attracting cycle. Such stable oscillations can be interpreted as a sign of overexploitation of the system (van Voorn et al., 2007). They are sometimes referred to as 'boomand-bust' cycles and have been observed episodically in ecosystems (Desta & Coppock, 2002); even leading to a successful prediction of the next crash in animal numbers (Coppock et al., 2008). As the management nullcline is shifted to the left, that is, toward V , the amplitude of the cycle gets bigger, until F I G U R E 4 Graphical summary of the general stability analysis, as a phase plane. Note that the vertical management nullcline g(V) − D ; (V) = 0 is not shown, as its position varies with . Decreasing (resp. increasing) moves the management nullcline to the left (resp. right). The phase plane is partitioned into five areas (labeled 1-5) where the vertical management nullcline can lie, and for each of which the system displays qualitatively different behavior. H reaches 0. We call "V " the critical value for V eq that corresponds to the stable cycle's destruction. Below V , stable limit cycles do not exist anymore. The destruction of this stable attractor corresponds to a global bifurcation, which marks the deterministic tipping from a productive system to the degraded state, that can be driven by changes in management.

regions (4) and (5)
These cases are undesirable, as the only stable steady state-and hence, the only attractor-is (0; 0). When the management nullcline is in these regions, the whole phase plane is the basin of attraction of the degraded state. This occurs when animals are not removed fast enough from the system or, put differently, the vegetation is too fragile for this type of exploitation.
Region (4) corresponds to the system collapsing due to overexploitation, whereas in region (5) the vegetation is unable to sustain itself even in the absence of herbivores. Bifurcations occur at B1, B2, and B3. Bifurcation B1 corresponds to the transition between regions 1 and 2 in Figure 4; bifurcation B2 corresponds to the transition between regions 2 and 3; bifurcation B3 corresponds to the transition between regions 3 and 4.

| Comparing the productivity of different destocking strategies
When considering livestock productivity, two main categories can be considered, namely systems of primary production, that involve slaughtering the animals for meat and other animal products, and systems of secondary production, in which the animals are kept alive while by-products such as dairy and wool are harvested. We consider the impact of varying the stocking density and adaptivity on these two production systems.
For a dairy production system with a given and any , we assume that the dairy production at equilibrium during the arbitrary time interval [0; t] is proportional to the number of animals so that production will be equal to where is a constant reflecting the per capita dairy production rate.
In the case of a meat production system with a given and any , the meat production at equilibrium during the arbitrary time interval 0; t is proportional to the sum of all the animals removed from the system during that interval, that is, assuming that all animals are slaughtered before their natural death.
We see that in both cases, the production output at equilibrium has no dependence on . Furthermore, in both cases, the production increases as the long-term stocking density H eq increases (i.e., as decreases).

| Comparing the resilience to drought of different destocking strategies
The resilience to degradation induced by management choices, that is, changes in the values of the parameters and , is measured by the metrics dist bif ( ) and dist bif ( ) (defined in Section 2), respectively. A first pathway to rangeland degradation occurs when the management choice of a long-term stocking density (dictated by the value of parameter ) leads to the deterministic collapse of the system. This corresponds to overexploitation-driven degradation: if the management nullcline is shifted to the left of V crit (in case of a subcritical Hopf bifurcation) or V (in case of a supercritical Hopf bifurcation), even in a gradual manner, then (0; 0) becomes the only stable steady state. Because the vegetation density V crit at which the Hopf bifurcation happens is independent of and , and because we know that decreasing shifts the equilibrium vegetation density V eq toward V crit , whereas does not affect V eq , we can conclude for the distance to bifurcation (for fixed and , respectively) that In other words, increasing the stocking goal H eq brings the system closer to a deterministic shift toward an oscillating or degraded system, whereas changing the adaptivity of the system has no effect.
The second theoretical pathway to rangeland degradation, which is always present, even in case of sensible management, is a consequence of perturbations to the state variables. When the parameters and initial conditions are such that the system is in a sustainable productive state V eq ; H eq , there is still a risk that external perturbations drive it to the stable degraded state (0; 0). This risk will vary, depending on the values of the management parameters and . The resilience of a system to an external perturbation consisting of a sudden loss of vegetation (as described in Section 2), as a function of the management parameters and , is measured through the resilience metric V eq V pmax ( ; ) (also described in Section 2). Deriving our results for this resilience metric takes the form of a mathematical proof. The interested reader can find explanations of unfamiliar concepts in any good linear algebra or introduction to dynamical systems textbook.
Readers who are uninterested in the mathematical proof can directly go from here to the summary of the results.
To derive the resilience metric V eq V pmax ( ; ), we need to study the threshold that separates the degraded and productive states' basins of attraction in the phase plane. In our case, this curve, the 'separatrix', is the stable manifold of the saddle point (A; 0), that is, On the other hand, because we are considering a sustainable productive system, we know that our management nullcline is in region 2 and, hence, to the right of the point (A; 0). Therefore, we have Now, we can derive results for our resilience metrics when is fixed and for h , l , such that h < l and such that their respective equilibria V eq h ; H eq h and V eq l ; H eq l are sustainably productive. We know from the previous subsection that H eq h > H eq l and V eq h < V eq l .

The stable eigenvectors are given by
where we have ignored the infinitesimally small parameter since A > 0 (by assumption). We have l V eq l > h V eq h , so the stable eigenvector is steeper in the system with l than in the system with h . Then,

Equation (3) becomes
As long as we are on the left-hand side of the management nullclines, the numerator is negative. Because we are considering trajectories coming into (A; 0), and because we know that below the hump-shaped vegetation nullcline the trajectories move away from (A; 0), all separatrices Γ ; necessarily lie above the vegetation nullcline. Hence, the denominator is also negative. We can then conclude that the ratio (4) is positive for i = h, l, and is greater for the system with the larger value. We have, therefore, shown for all V on the left-hand side of the management nullcline, that is, the separatrix of the system with l is always on the "outside" of the separatrix of the system with h on the left-hand side of the management nullcline. We can, hence, conclude, for any fixed and any h , l producing stable productive equilibria and such that h < l , that that is, the threshold for vegetation density before collapse is lower for a greater , therefore that is, a productive system with greater can recover from greater percentages of vegetation loss.

Now that we have compared the resilience of a sustainable
productive system for a fixed and different values of , we consider the case where is arbitrarily fixed such that it yields a sustainable productive system, and we have c and a such that a > c . Noticeably, following our definition, for a fixed , all values of yield the same management nullcline and hence the same longterm stocking goal H eq . Following the same reasoning as previously, we find that that is, the threshold for vegetation density before collapse is lower for a greater , therefore that is, a productive system with greater can recover from greater percentages of vegetation loss.

| Summary of the results for the maximal loss ratio
We have shown that, for any fixed and any h , l producing stable productive equilibria and such that h < l , that is, the threshold for vegetation density before collapse is lower for a greater , therefore that is, a productive system with greater can recover from greater percentages of vegetation loss.
We have also shown that, for any fixed producing a stable productive equilibrium and any c , a such that c < a , that is, the threshold for vegetation density before collapse is lower for a greater , therefore that is, a productive system with greater can recover from greater percentages of vegetation loss.
In conclusion, on the one hand, a greater value (and therefore, a lower long-term stocking goal H eq ) provides greater resilience to the sustainable productive state, when considering perturbations in the negative-V direction. We can conclude that for any given degree of adaptivity , systems with lower long-term stocking goal H eq have greater resilience in virtue of the expected trade-off between resilience and productivity in grazing systems. On the other hand, we have also demonstrated that greater adaptivity of the destocking rate provides greater resilience of the productive system to vegetation losses. These two results mean that increasing or will systematically decrease the risk of collapse to a degraded state, as illustrated in Figures 6 and 7.

| Illustration of the results
We generated Figures 6 and 7 using the realistic model specification and parametrisation given in the appendix. This example is derived from the well studied Klausmeier model (Klausmeier, 1999) for dryland vegetation, which consists of a system of two partial differential equations, one for plants and one for water, as explained in the appendix.
As noted earlier, the choice of parameters is for illu stration purposes only and has no bearing on the general results shown above, which are valid for all parameters.
However, we note that the value of H max yielded by our model, We compare numerically the behavior and resilience of this specific system, using the pplane8 phase plane plotter package of Matlab, (Harvey, 2022), with the Dormand-Prince (ode45) solver, for: • two different values (at fixed = 0): we pick h and l such that they, respectively, yield the stocking targets H eq h = 0.9H max and H eq l = 0.6H max (Figure 6).
• two different values (at fixed = h , as defined just above): we pick c = 0 and a = 1, yielding, respectively, what we call "constant" and "adaptive" strategies ( Figure 7).
We apply to all four scenarios two perturbations with different magnitudes, namely the weaker perturbation pert 1 and the stronger perturbation pert 2 . In Figure 7, pert 1 is defined as the loss of 50% of the vegetation at equilibrium, whereas pert 2 is defined as the loss of 66% of the vegetation at equilibrium. In Figure 6, the equilibria for l and h do not coincide, and pert 1 (respectively, pert 2 ) is defined as the loss of 50% (respectively, 66%) of V eq h , the vegetation density at equilibrium associated to h . Therefore, in the higher case, l , these perturbations are even more severe, which strengthens our results according to which higher values provide more resilience to sudden losses of vegetation.
Each subfigure is the superposition of the phase planes of the two strategies we are comparing. In all cases, the vegetation nullcline ℋ(V) (dotted pink) has the same position, since it is affected by neither the value of nor α. When is fixed, the management nullclines (dotted yellow lines) of the two parametrisations also coincide (Figure 7). The full green dots represent the stable steady states. Again, they coincide if and only if the parametrisations feature the same value for . The dashed lines are the separatrices, that is, the limits of the basins of attraction of the sustainable productive states. The full arrowed lines represent the trajectories of the system after a perturbation. In Figure 6, we are comparing the grazing systems' trajectories for h and l ( = 0 is fixed). We note that the greater value l yields a lower stocking target H eq l and that the basin of attraction for l englobes the one for h . In subfigures (a) and (b), the equilibrium points V eq l ; H eq l and V eq h ; H eq h are marked with a green circle containing the letters 'l' and 'h', respectively. In subfigures (a), (c), and (e), we see that after a perturbation consisting of halving the V eq h vegetation density (pert 1 ), both systems recover to their respective equilibria. In both cases the system was still inside the basin of attraction of its productive state after the perturbation, as can be seen in subfigure (a). In subfigures (b), (d), and (f), after a perturbation consisting of dividing by three the V eq h vegetation density (pert 2 ), only the system with greater destocking parameter l and lower stocking target H eq l recovers. The system with h ends up degraded as the perturbation kicked the system out of the basin of attraction of V eq h ; H eq h . We see in (b) that the shape of the two separatrices in the direction of perturbation does not differ widely and we can infer from the plots (d) and (f) that the system with l owes its recovery to its lower stocking target H eq l .
This illustrates the trade-off between resilience and productivity in grazing systems.
In Figure 7, we are comparing the constant ( c = 0) and adaptive ( a = 1) destocking strategies, for a fixed = h . In both cases the system has the same stable productive equilibrium V eq ; H eq and, therefore, the same stocking target and productivity. In subfigures (a), (c), and (e), after a perturbation consisting of halving the V eq vegetation density, both systems recover to V eq ; H eq . In both cases the system was still inside the basin of attraction of its productive state after the perturbation. In the case of the adaptive strategy a , we can tell from the shape of the trajectory in (e) that the herd size has fluctuated more, due to a faster destocking rate after the perturbation.
This could appear as a downside of this strategy. However, in sub- , and (f), after a perturbation consisting of dividing by three the V eq vegetation density, only the adaptive system recovers.

F I G U R E 6
Comparing the resilience of a sustainable productive system to the effects of drought for different baseline destocking rates l and h , where h < l and with fixed adaptivity = 0. Left panels (a, c, e) show the response to a weak perturbation while right panels (b, d, f) show the response to a stronger perturbation. Top panels (a,b) show the phase-space while the middle and bottom (c, e, d, f) show the trajectories of the biomass variables over time.  which is a proxy for productivity. The variation in productivity shown on the horizontal axis corresponds to a variation in : greater baseline destocking rates lead to lower long-term stocking densities H eq and conversely. The negative slope of all three curves depicts how F I G U R E 7 Comparing the resilience of a sustainable productive system to the effects of drought for different degrees of adaptivity c = 0 and a = 1, where the baseline destocking rate = h and hence the long-term stocking density H eq is fixed. Left panels (a, c, e) show the response to a weak perturbation while right panels (b, d, f) show the response to a stronger perturbation. Top panels (a,b) show the phase-space while the middle and bottom (c, e, d, f) show the trajectories of the biomass variables over time. imal when the productivity is at herbivore carrying capacity H max , that is, at the system's Hopf bifurcation. We saw in the general results (Section 3) that the resilience decreased as H eq increased, which can be observed for all three curves. We also saw that the resilience increased as adaptivity increased, which is illustrated through the improvement in resilience as is increased.

| DISCUSS ION
A main result from this work was that dry rangelands with coupled plant-herbivore dynamics can benefit greatly, in terms of their resilience and without loss of long-term productivity, from adaptivity of the stocking density to the lack of forage. We found that both a low stocking density and a high stocking adaptivity, achieved through modulation of the destocking rate, provided crucial resilience in dry rangeland systems subject to drought.
Our modelling approach was very general and relied on simple and realistic constraints. The analysis is, therefore, valid for a variety of more specific models that fall under its wide umbrella. The results are robust and independent from parametrisation and from specific and possibly arbitrary choices of functional forms to represent specific mechanisms.
For any given strategy, our model results showed that a higher long-term stocking target implies a lower resilience to sudden vegetation losses. This trade-off between productivity and resilience is well known and observed in various exploitation models (Fletcher & Hilbert, 2007;Noy-Meir, 1975) and is intuitively justified by the additional stress put on vegetation by a higher number of animals.
Mathematically, in our work, this translated into the basin of attraction of a productive state contracting as the long-term stocking target increases. Consistent with earlier work on general herbivorevegetation systems (Noy-Meir, 1975), we found that there exists a maximal sustainable livestock density, that is, an animal carrying capacity, H max . The existence of an animal carrying capacity and of a trade-off between productivity and resilience is consistent with recommendations for low stocking densities, as stocking well below carrying capacity lessens the risks of degradation following a drought (McLeod, 1997;Vetter, 2005). This trade-off also highlights the potentially deleterious effects on resilience of strategies aiming to increase productivity, such as feed supplementation (Müller et al., 2015), or improvement in veterinary care, which decreases animal mortality in times of drought and, therefore, reduces the degree of coupling between vegetation and animals. A better understanding of these effects is essential to prevent counterproductive policies or practices. Maximizing resilience instead of short-term productivity is also in accordance with recommendations from current frameworks of food system resilience (Rotz & Fraser, 2015).
Policy makers for sustainable rangeland management have often favored fixed stocking densities rather than traditional adaptive pastoralism, which has been perceived as a source of degradation (Lv et al., 2019). Our destocking strategies with > 0 fell into the category of adaptive stocking strategies, as they regulated the number of animals in reaction to perturbations in the vegetation density. We found two results supporting the use of adaptive strategies in dry rangelands. First, we found that the animal carrying capacity H max did not vary with different degrees of adaptivity of the destocking rate: no matter how sensitive the destocking rate was to vegetation changes, all strategies allowed the same maximal number of animals in the system. Second and importantly, a higher dependence of the destocking rate on the lack of forage systematically increased the resilience of the productive state, while maintaining the long-term production goal. This means that adaptivity in the destocking rate enabled circumventing the usual trade-off between productivity and resilience. This result is in line with previous, less general, studies highlighting the importance of adaptivity in dry rangelands (Freier et al., 2014;Jakoby et al., 2015). As suggested by van de Koppel and Rietkerk (2000) and unlike Jakoby et al. (2015), we did not find that adaptivity eliminated the risk of a catastrophic collapse due to overstocking. We found that the higher the adaptivity, the closer resilience was to the no-herbivore system's resilience, which itself is susceptible-by construction of our model-to irreversible droughtdriven collapses.

F I G U R E 8
Relation between the resilience and the productivity of a sustainable productive system, for different values of the management parameters and . On the horizontal axis, the long-term stocking density H eq ( ) is a proxy for the rangeland's productivity. On the vertical axis, the resilience is measured through the maximal loss ratio, Our results also contributed to the theoretical debate on dry rangeland modeling, which opposes supporters of the equilibrium and non-equilibrium theories. The latter argue that even though a carrying capacity exists, the frequency and intensity of droughts routinely force the system away from it and, therefore, make it irrelevant. In this framework, dry rangelands are, therefore, considered not to be vulnerable to overgrazing-driven degradation, precisely because the animal density will drop after a drought (Ellis & Swif, 1988), due to decreased vegetation but also increased thirst and illness related mortality (Catley et al., 2014). This argument agrees thoroughly with our result: increased adaptivity, that is, a sharp decrease in the number of animals following a drought, will make dry rangelands more resilient to drought.
Our general results apply to dry rangelands, but they can also be applied to other consumer-resource (or even predator-prey) systems, if the consumer's (or predator's) loss rate depends on the amount of resources available. The realistic assumption that resource shortage is correlated to consumer's weakening and lower resistance to diseases, which lead to extra mortality (Catley et al., 2014), is not accounted for in other models where, instead, the consumer loss rate is constant. We are aware of only one other predator-prey model where the predator mortality rate depends on prey density (Minter et al., 2011). In that work, which concerns protozoa, a prey density dependent death rate is derived experimentally for the predator and leads to quantitatively significant differences in the resulting dynamics. The authors show the importance of taking into account such variability in the predator's (consumer's) death rate to improve population models. In our case, we saw that resource-dependent consumer loss rate can be a management tool that enhances the system's resilience.
Our study contained several limitations and shortcomings that could be addressed in future research. First, our model was a meanfield representation of a system that is spatially extended in the real world. Not modeling explicitly the spatial dimension prevented our system from displaying spatial organization characteristics that are crucial in many dryland ecosystems. For instance, spatially explicit mathematical models of dryland vegetation elucidate how dryland vegetation, characterized by local plant-plant facilitation mechanisms and longer-range competition for water, self-organizes into patterns (Kéfi et al., 2007;Klausmeier, 1999). Importantly, according to theoretical work by Rietkerk et al. (2021), this patterning allows vegetation to evade the abrupt degradation predicted by our meanfield model, by undergoing a more gradual change instead. However, when modeling herbivory, it has been shown that grazing can attenuate the buffering provided by spatial patterning against sudden degradation (Siero et al., 2019). Therefore, it is unclear how incorporating the spatial organization of dryland vegetation would affect our results.
Accounting for the spatial dimension would also allow for the representation of animal mobility. Nevertheless, our model already enabled us to infer that mobility was beneficial in terms of resilience to drought: resilience increased if animals disappeared from the system when vegetation density was low. This disappearance could correspond to animal displacement. The animal density could then be restored once vegetation had recovered, which would correspond to the removed animals coming back or new animals being restocked.
We could adapt our current model to better represent pastoral mobility and the practices of transhumance or rotational grazing by explicitly modeling space in a discrete manner, by defining several different pastures. Such an approach would take into account the dynamics of the displaced animals and of the alternative pasture(s), whose coupling could affect the dynamics of the original pasture, for example by acting as key resources (Illius & O'Connor, 1999;von Wehrden et al., 2012). A study on rotational grazing (Chen & Shi, 2018) finds that both production yields and stockpiled forage increase for many rotational configurations. Interestingly, there is an increasingly positive perception of pastoral mobility among dryland researchers (Adriansen, 2005) and recent studies aim at incorporating traditional indigenous knowledge into policy making (Selemani, 2020). This is consistent with recent modeling research showing that grazing, when managed in a spatially non-uniform manner, can improve resilience to droughts (Zelnik et al., 2021).
Future research might also incorporate different modalities of environmental perturbations. Our model used a constant rainfall rate, which allowed us to focus on fluctuations in herd size that were driven by animal-plant feedback rather than by rainfall stochasticity. However, this assumption is unrealistic, as real-life (semi-)arid regions usually experience high variability in their rainfall rates, and these play quite some role in their vegetation dynamics (Baudena et al., 2007;Verwijmeren et al., 2021). Still, even though we did not directly model intra-and interannual variation in rainfall and hence in vegetation growth, the mathematical analysis of our model proved that if random pulse perturbations are applied, then lower stocking densities and higher adaptivity of the destocking rate systematically decrease the probability of degradation. Further work could study sequences of droughts that allow only partial recovery between the perturbations, or press perturbations, where the duration of the perturbation is prolonged in time, rather than modeled as an instantaneous event.
Finally, our study lacked an economic dimension since it did not take into account destocking costs or market considerations. Even though adaptive destocking did not imply any long-term reduction in productivity, we can expect that in reality it would be associated with increased logistic costs. A more realistic bioeconomic model would need to consider market fluctuations as well as the transient production linked to extra destocking, when computing the productivity of a system. For example, a case study in the Sahel reveals how livestock owners' decision to sell-and hence their ecological impactdepends on their type of access to the market, as well as complex institutional and cultural factors Turner and Williams (2002). This study emphasizes the need to understand and incorporate local livestock market specificity when designing rangeland management policies in drought-prone areas.
Water-scarce and drought-prone rangelands are vulnerable ecosystems that were traditionally managed sustainably thanks to livestock mobility (Freier et al., 2014). Now, in times of increased sedentarisation and intensification of stocking densities (Reid et al., 2014), understanding the positive impact of adaptivity on the resilience of vulnerable systems can help in the design of sustainable management strategies and policies. In practice, our study supports initiatives that facilitate adaptive destocking actions and/ or livestock mobility. Initiatives to facilitate adaptive destocking actions already exist and include facilitating slaughtering at the onset of drought, easing access to the market, and transforming and giving value to destocked meat, for example, implementing logistics to dry and/or can destocked meat and distribute it as a supplement in times of drought (Abebe et al., 2008;Morton & Barton, 2002). Sherratt: Conceptualization (supporting); formal analysis (supporting); methodology (supporting); supervision (equal); writing -review and editing (equal).

ACK N OWLED G M ENTS
We are very grateful to the anonymous reviewer for their outstanding contribution. They have provided insightful comments that have significantly improved the paper and spotted a mistake present in the original manuscript. We would also like to thank

UK Research and Innovation-Engineering and Physical Sciences
Research Council (UKRI EPSRC) grant EP/S023291/1; Heriot-Watt University; University of Edinburgh.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.

Model example
To illustrate our general model, we provide a specific example that satisfies the general form and all constraints on (1): where V is the aboveground biomass density, j the yield of plant biomass per unit water consumed, r the rate of water uptake by plants, m the plant mortality rate, R the constant rainfall rate and e the rate of evaporation and loss to the deeper ground, the maximal consumption rate per herbivore and the half persistence parameter.
It explicitly captures how plants locally increase water availability compared with bare ground and features implicit competition for water. The non-zero roots of f(V) are which illustrate the dependence of the vegetation carrying capacity K and Allee threshold A on environmental and plant-specific parameters. In particular, here A and K are, respectively, a decreasing and an increasing function of the average rainfall rate R.

A.1 | Derivation of the model
The spatially uniform version of Klausmeier model equation (Klausmeier, 1999) is Assuming that root-water infiltration dynamics occur on a much faster time-scale than vegetation growth/decay and grazing dynamics, we can apply a quasi steady state approximation. Hence we suppose that available water is always at equilibrium so which yields

Substituting this in the vegetation equation yields
By adding coupled herbivores dynamics where we define g(V) as the widely used Holling II functional response, we have the resulting example system (5).
All parameters are fixed, taken or derived from the literature, as given in Table A1, with the exception of and , our control parameters. We are considering smallstock, that is, sheep or goat ("shoat") animal units, such that the whole flock is female (negligible number of males) and each female lambs one offspring per year on average when vegetation is at carrying capacity (V = K). Hence, the maximal per capita growth rate is 2 year −1 , and we can solve for .
The mean annual precipitation rate R is chosen to be relatively low but still above the threshold for equilibrium dynamics in our coupled system of ODEs (Boone & Wang, 2007;Briske et al., 2003).